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ABSTRACT 

The recent ALMA observations of the disk surrounding HL Tan reveal a very complex dust spatial 
distribution. We present a radiative transfer model accounting for the observed gaps and bright rings 
as well as radial changes of the emissivity index. We find that the dust density is depleted by at least a 
factor of 10 in the main gaps compared to the surrounding rings. Ring masses range from 10-100 M 0 
in dust, and we find that each of the deepest gaps is consistent with the removal of up to 40 M 0 of 
dust. If this material has accumulated into rocky bodies, these would be close to the point of runaway 
gas accretion. Our model indicates that the outermost ring is depleted in millimeter grains compared 
to the central rings. This suggests faster grain growth in the central regions and/or radial migration 
of the larger grains. The morphology of the gaps observed by ALMA - well separated and showing 
a high degree of contrast with the bright rings over all azimuths - indicates that the millimeter dust 
disk is geometrically thin (scale height ~ 1 AU at 100 AU) and that a large amount of settling of large 
grains has already occurred. Assuming a standard dust settling model, we find that the observations 
are consistent with a turbulent viscosity coefficient of a few 10“^. We estimate the gas/dust ratio 
in this thin layer to be of the order of 5 if the initial ratio is 100. The HCO+ and CO emission is 
consistent with gas in Keplerian motion around a 1.7 Mq star at radii from < 10 — 120 AU. 

Subject headings: stars: individual (HL Tan) - protoplanetary disks - stars: formation - submillimeter: 
planetary systems - techniques: interferometric - radiative transfer 


1. INTRODUCTION 

Planet formation within disks around young stars re¬ 
quires the growth of sub-micron dust grains over many 
orders of magnitude. The dust feeding the disk in its 
early evolution is thought to have grown from the sub¬ 
micron grains typical of the interstellar medium, up to 
micron sized particles in the dense regions of molecular 
clouds and cores (e.g. Pagani et al. 2010). Once incor¬ 
porated into the disk, these grains must then grow from 
micron-sized to pebbles and ultimately kilometer-sized 
bodies. This extensive growth phase can only occur in 
the highest density regions of a protoplanetary disk - at 
or close to the disk’s midplane (Dominik et al. 2007; Testi 
et al. 2014). 

The gas drag acting on the dust as it orbits the cen¬ 
tral star causes it to settle rapidly to the midplane. Just 
how much remains a poorly constrained quantity by di¬ 
rect observations. The small dust particles (< 1/im), 
best traced by scattered light, remain coupled tightly to 
the gas, while the larger particles, best studied at longer 
wavelengths, move toward the midplane. The thickness 
of the midplane depends on the balance between verti¬ 
cal settling and counter-acting stirring mechanisms (Ga- 
raud et al. 2004; Dullemond & Dominik 2005; Fromang 
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& Papaloizou 2006). Once in the midplane, particles 
are also expected to drift radially inward as a result 
of the combined action of the stellar gravity, the gas 
drag and the radial pressure gradient (Barriere-Fouchet 
et al. 2005; Birnstiel et al. 2010). The presence of plan¬ 
ets within the disk further complicates the dynamics of 
the dust grains as they may create gaps whose profile 
depends on the planet mass, density profile, and grain 
size (Paardekooper & Mellema 2004; Fouchet et al. 2007; 
Gonzalez et al. 2012). 

Fitting multi-color scattered light images of edge-on 
disks indicates the presence of stratification, with smaller 
grains found higher up in the disk atmosphere while 
larger grains have settled deeper (e.g. Pinte et al. 2007, 
2008b; Duchene et al. 2010), as predicted by models. 
Tracing the properties of the dust layers deep inside 
the disk requires observations at long (millimeter) wave¬ 
lengths where the disk becomes optically thin. Several 
multi-wavelength studies are now suggesting that the 
larger grain abundances depend on radius, with a concen¬ 
tration toward the center of protoplanetary disks (Ricci 
et al. 2010; Guilloteau et al. 2011; Perez et al. 2012; Menu 
et al. 2014); however, the spatial resolution was insuffi¬ 
cient to resolve the vertical structure of the dust layer. 
In general, the midplane region of disks remains poorly 
constrained by direct observations. 

HL Tan is a T Tauri star of spectral energy distribu¬ 
tion (SED) Class I-H in the Taurus molecular cloud, at 
a distance of 140 pc. The disk has been widely studied 
in the (sub-)millimeter wavelength continuum (Beckwith 
et al. 1990; Mundy et al. 1996; Wilner et al. 1996; Lay 
et al. 1997; Chandler & Richer 2000; Looney et al. 2000; 
Greaves et al. 2008; Kwon et al. 2011; Stephens et al. 
2014). The highest resolution millimeter observations 


2 


so far (Kwon et al. 2011) revealed a 120 AU dust disk 
and suggested, with comparison with the SED, that the 
millimeter grains may have already settled to the mid¬ 
plane. In addition to the disk structure, an orthogonal 
optical jet and a molecular bipolar outflow (e.g. Mundt 
et al. 1990) and an envelope (Men’shchikov et al. 1999; 
Robitaille et al. 2007) have been reported. 

HL Tau was recently observed with ALMA using base¬ 
lines up to ^ 15 km at wavelengths from 0.8-3 mm, which 
resulted in spatial resolutions down to 3.5 AU (ALMA 
Partnership et al. 2015b, hereafter ALMA15). These 
data provide the critical angular resolution needed to 
study the details of the disk midplane. They indicate 
the presence of numerous gaps and rings. 

Here we model the recent ALMA data at all three 
wavelengths, in order to reproduce the main features of 
the HL Tau disk. We interpret the results in terms of the 
underlying dust distribution and evolution. We investi¬ 
gate whether a single unified dust model can reproduce 
both the ALMA images and SED, and compare the dust 
and gas disk also observed by ALMA. 

2. OBSERVATIONS 

HL Tau was observed as a science verification target 
within the ALMA Long Baseline Campaign; full details 
of the observations are given in ALMA15 and ALMA 
Partnership et al. (2015a). The band 3, 6, and 7 cali¬ 
brated continuum data sets and CLEANed images were 
downloaded from the ALMA Science portal^, together 
with the tables for self-calibration, which we applied to 
the data. 

Erom the images, the observed brightness in each band, 
and the spectral index profiles measured along the disk 
major axis are presented in Eigure 1. A similar struc¬ 
ture of alternating bright rings and dark gaps is seen in 
all three bands. The spectral index measured between 
0.87 and 2.9 mm (lower panel) displays a significant in¬ 
crease from the inner regions to the outer disk, from ^2.2 
to ^3.5, indicating a radial change in the spatially con¬ 
volved dust emissivity and/or optical depth. The com¬ 
bined cleaning of the band 6 and 7 data (ALMA15) yields 
a higher resolution spectral index map, and a cut through 
this along the major axis is presented in the middle panel 
in Eigure 1. As noted previously, this shows that the op¬ 
tical depths are significant in the rings at these shorter 
wavelengths, with sharp increases in the spectral index 
inside the gaps. In order to extract the underlying dust 
structure, we reproduce the data at all three wavelengths, 
and then compare the models with these radial profiles. 
Any differences in the derived models would then suggest 
evidence of a radial change in the dust characteristics. 

The calibrated spectral line data sets in band 3 for 
J=l-0 and HCO+ J=l-0 were also downloaded 
and CLEANed. The default CLEAN parameters and uv 
taper from ALMA15 were used for HCO+ giving a res¬ 
olution of 0.25arcsec. Eor the CO data set, CLEANing 
with a uv taper giving a spatial resolution of 0.2arcsec 
was found to be a good compromise between resolving 
the disk and obtaining enough signal/noise per beam. 
Although the SNR of the spectral line data sets is rel¬ 
atively low and suffers from contamination at ambient 
cloud velocities, we were able to compare the results with 

^ http://www.almascience.org/ 



Figure 1. Top: cuts of observed and synthetic CLEAN maps 
along the major axis of the disk (red = band 3, 2.9 mm, green = 
band 6, 1.3 mm, blue = band 7, 870p,m.). The spatial resolution 
along the disk is 10, 4, and 3AU in the three bands. Middle: 
spectral index profile derived from the combined 6+7 bands spec¬ 
tral index map, from ALMA15. Bottom: Spectral index profile 
obtained by convolving the band 6 and 7 maps to the same resolu¬ 
tion as the band 3 map and fitting a spectral index pixel-by-pixel. 
In all three panels, observations are plotted with a full line while 
the model is plotted with a dashed line. 

a simple Keplerian disk in order to investigate the gas 
distribution. 

3. CONTINUUM MODELING APPROACH 

To model the ALMA continuum data, we use the radia¬ 
tive transfer code MCEOST (Piute et al. 2006, 2009). Be¬ 
cause of the remarkable symmetry of the ALMA images, 
we assume an axisymmetric, i.e. two-dimensional (2D), 
density structure. In particular, we ignore the slight off¬ 
sets and change of inclinations of the rings compared to 
the central star (ALMA15). We use standard disk pa¬ 
rameters to describe the disk structure {e.g. Pinte et al. 
2008a; Andrews et al. 2009, see Section 3.1) except for 
the surface density which shape is let free and is fitted 
for. Our aim is to determine whether such a simple model 
of the disk density structure can reproduce the main fea¬ 
tures observed by ALMA. Due to the complexity of the 
ALMA data sets, we build the model of surface density 
from the images and checked a posteriori by calculating a 
that the model and observed visibilities were in good 
agreement. 

3.1. Radiative transfer model 

We assume a flared gas density structure with a Gaus¬ 
sian vertical profile pg{r,z) = pQ{r) ex.-p{—z^/2 h{r)^) 
where r is the distance from the star in the disk mid¬ 
plane and 2 ; is the altitude above the midplane. Eor the 
gas disk scale height h(r), we use a power-law distribu¬ 
tion h{r) = ho (r/ro)^ where ho is the scale height at 
radius ro = 100 AU. To keep the problem tractable, we 
fix some parameters and explore the parameter space 
in terms of surface density and dust settling. The disk 
scale height is set to 10 AU to be consistent with models 
of thermally supported disks {e.g. Close et al. 1997). P 
is the disk flaring exponent, which we set to 1.15, again 
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consistent with typical disk models. The exact values 
of ho and P do not strongly affect the final conclusions 
and the defined gas disk structure is very close to hydro¬ 
static equilibrium (see Appendix B). The inner radius is 
set to the dust sublimation radius. The outer radius is 
set to 150 AU, i.e. large enough that the surface density 
reaches negligible values at this distance from the star 
and in agreement with the size of the observed disk. 

A surface density distribution per grain size i;(r, a) is 
then fitted at each wavelength independently, using the 
procedure described below (Section 4.1). 

Dust settling is implemented by describing the dust 
transport by turbulence as a diffusive process. This was 
shown to be a reasonable approximation to global MHD 
simulations of stratified and turbulent disks (Fromang 
& Nelson 2009). The degree of dust settling is set by 
varying the turbulent viscosity coefficient a (Shakura & 
Sunyaev 1973). We assume that the gas vertical profile 
remains Gaussian and that the diffusion is constant ver¬ 
tically and given hy D = Cgh a/ Sc^ where = h D is the 
midplane sound speed, h is the gas scale height, and Sc 
is the Schmidt number, which we fix to 1.5. The vertical 
density profile for a grain of size a is given by equation 
19 of Fromang & Nelson (2009): 


p(r, z, a) (X S(r) exp 


firs(a) f ^ 

e2h,2 1 


D V 

J 2h?\ 


( 1 ) 


with D = a /D = and rs(a) = the dust 

stopping time, where pd is the dust material density and 
pg is the gas density in the midplane. The density is 
finally normalized so that 



dz = i;(r) 


( 2 ) 


We assume passive heating and fix the properties of 
the central star to Teff = 4 000K and = 11 Tq 
(M en’shchikov et al. 1999). We use a central mass 
= 1.7 Mq, as derived from the CO and HCO+ obser¬ 
vations (see Section 5). 

We adopt a fixed dust mixture composed of 60 % sil¬ 
icate (Dorschner et al. 1995), 15% amorphous carbon 
(Zubko et al. 1996), and 25% porosity, although the de¬ 
tails of the dust composition do not affect the results 
presented in this paper. We use a grain size distribution 
dn(a) oc aPda between amin and Umax- <^min and Umax are 
set to 0.03/rm and 3 mm and we initially use p = —3.5 
(integrated over the whole disk). However, note that 
Umax and p are both modified locally due to the dust set¬ 
tling and varying surface density as a function of grain 
size. Each grain size is represented by a distribution of 
hollow spheres with a maximum void fraction of 0.8. 

In summary, we fixed most of the model parameters, 
and focused our analysis on the parameters that are 
mainly probed by the new ALMA data: the disk sur¬ 
face density per grain size E(r, a) and dust settling via 
the a parameter. 


3.2. Building a model of the HL Tau disk 

The construction of the 2D disk density structure was 
performed in two successive steps. First, we build the 


surface density of the model by reproducing the bright¬ 
ness distribution along the major axis, where the line- 
of-sight confusion due to the projection of the vertical 
structure is minimal. This procedure is repeated for the 
three ALMA bands and the resulting surface densities, 
which correspond to various grain sizes, are combined in 
a single model. In a second step, the vertical distribu¬ 
tion is explored by modeling the complete 2D data sets 
and quantitatively comparing the synthetic and observed 
visibilities. The final model shows residuals smaller than 
the map rms, and visibilities in very good agreement, in¬ 
dicating that the model resulting from this procedure can 
account for most of the observed features in the ALMA 
data sets. 

To obtain a surface density profile that matches the 
observed images, we use an iterative procedure. We first 
extract a brightness profile along the major axis of the 
disk in the CLEANed map and take an average of both 
sides. Because the model will be convolved before com¬ 
parison with observations, we deconvolve the extracted 
brightness profile by a Gaussian beam whose EWHM cor¬ 
responds to the projected EWHM along the disk semi 
major-axis of the observed 2D beam. This deconvolved 
brightness profile is then converted to a surface density 
profile assuming an initial power-law radial temperature 
profile (exponent = -0.5) consistent with the radial tem¬ 
perature dependence in typical disk models (D’Alessio 
et al. 1998; Dullemond et al. 2002; Pinte & Laibe 2014). 
The assumption of power-law temperature with expo¬ 
nent -0.5 is needed only for to the first iteration, T(r, z) 
is calculated for all subsequent iterations. Using this 
initial surface density profile, we compute a model us¬ 
ing MGEOST, convolve it with the corresponding beam, 
and extract a model brightness profile along the major 
axis. The surface density at any point in the disk is then 
corrected by the ratio of the observed profile divided by 
the model profile. To avoid too many oscillations in this 
procedure, this ratio is limited between 0.8 and 1.2. The 
procedure is iterated until the model brightness profile 
does not change by more than 1 % at any radius. The 
procedure requires of the order of 30 iterations to reach 
convergence. 

Note that the procedure breaks down in the very inner 
regions of the disk due to high optical depth and be¬ 
cause the central beam encloses all the azimuthal angles. 
Eor the first five central AU, we fixed instead a power- 
law surface density profile. We find that a slope of -0.5 
reproduces well the brightness distribution of the outer 
parts of the central peak, but the exact value is not well 
constrained because the region is only partially resolved. 

Despite exploring a range of parameter space for both 
the disk geometry and the dust properties, we could not 
reproduce the compact and peaked emission at the po¬ 
sition of the star. Because the emission in the central 
peak is optically thick, we conclude that an additional 
unresolved source such as coronal free-free emission or 
extra heating due to accretion might be present close to 
the the star. We added the required unresolved emission 
at the position of the star in the model maps: 1.57, 1.90, 
and 2.59 mJy at band 3, 6 and 7 respectively. This rep¬ 
resents 2%, 0.2% and 0.1% of the total flux in bands 
3, 6, and 7 respectively. Note that the relative contribu¬ 
tion of this additional emission in the central beam is not 
well constrained and depends on the model parameters: 
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Figure 2. Upper: comparison of the observed and model clean maps in band 6+7. The scale units are mjy.beam”^. Note the different 
scale for the residual panel, which ranges from — Icr to +1(7, where a is the standard deviation in the combined image. Lower left: 
comparison between the observed (full lines) and model (dashed lines) visibility amplitudes along the disk major axis in the first spectral 
window of each band (band 3: red, band 6: green, band 7: blue). The model visibilities have been obtained from the MS file generated 
by the CASA simulator using the exact same antenna positions as in the observed data. The bins were constructed to all have the same 
number of visibility points. All the uv data points within 30° of the position of the major axis where considered. Lower right: comparison 
between the observed (red) and model (blue) real parts of the visibilities in the first spectral window of band 7 (337.494 GHz; other spectral 
windows give similar results). The inset is a zoom on the baselines larger than 500kA. The full lines correspond to the disk minor axis 
while the dashed lines correspond to the major axis (uv points within 30° of each axis were considered). 


the non-resolved contributions cannot be large, but their 
exact value should be taken with caution. 

For the iterative procedure, model images were pro¬ 
duced by convolving the MCFOST images by a Gaussian 
beam of equivalent FWHM. The final simulated visibil¬ 
ities and images were produced by applying a modified 
version of the CASA simulator to the model. This com¬ 
putes the model visibilities at the exact same (u^v) coor¬ 
dinates as the observations using the ft task. To compute 
the model visibilities, we first average the observed data 
in time (using 300s bins) and in frequency to obtain one 
single continuum channel per available spectral window. 
This was performed to reduce the number of (u^v) points 
and corresponding computing time of the CASA simula¬ 
tor. The thermal noise at each wavelength was set using 
the median value of the precipitable water vapor (PWV) 
of the observations: 1.3, 0.65, and 0.54 mm at bands 3, 6, 
and 7 respectively. The resulting data set of model visi¬ 
bilities was CLEANed using the same parameters as for 
the observations, and the observed and modeled images 
were compared for consistency (Figure 2). We estimate 
the quality of the model by calculating a between 
the observed and synthetic visibilities. After binning, we 
have 229 910, 420 668 and 571960 visibilities in bands 


3, 6, and 7 respectively, were computed using the 
formal uncertainties from the interferometric weights, to 
which we added quadratically a 10 % calibration uncer¬ 
tainty. The from the three individual bands bands 
were added to compute the final x2. Our final model has 
a total of 1282. For the final model, visibilities were 
also computed without performing any time averaging, 
and we checked that the reconstructed simulated images 
almost did not show any difference, and that the were 
not affected by more than 2 %. 

4. RESULTS FROM MODELING 

4.1. Surface density profile and radial concentration of 
large grains 

Figure 3 shows the resulting surface densities obtained 
by inverting the brightness profiles of the best models 
for the different wavelengths. The surface density pro¬ 
files derived from observations in the three ALMA bands 
are broadly consistent with each other, suggesting that 
our basic modeling assumptions are also correct. In par¬ 
ticular, we find that the central peak as well as the two 
most central rings are marginally optically thick, even at 
2.9 mm (Figure 3), while the outermost broad ring be¬ 
comes optically thin. The surface density extracted from 
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the Band 6 data is also in good agreement with the ana¬ 
lytical surface density derived by Kwon et al. ( 2011 ) and 
Kwon et al. (2015) from CARMA observations at the 
same wavelength (green dashed line in Figure 3, in the 
bright rings, the gaps were not detected by CARMA). 
In particular, our procedure recovers the transition be¬ 
tween the power-law description and exponential taper¬ 
ing at radii larger than 80 AU, as well as the shape of 
surface density in both regimes. 

Interestingly, assuming a grain size distribution slope 
of p = —3.5, the surface density profiles we derive from 
the three bands start to differ significantly in the outer¬ 
most ring, where the estimated surface density decreases 
as a function of wavelength. A model with a steeper 
slope, p = —4.5, results in similar surface density profiles 
at all bands in the outer ring, but the densities now dif¬ 
fer significantly in the central regions. This suggests that 
grains larger than a fraction of a millimeter are depleted 
in the outer ring compared with the rest of the disk, sup¬ 
porting the idea that grain growth is less efficient in this 
part of the disk or that in the event of grain growth, 
radial migration has already effectively removed a signif¬ 
icant fraction of the large grains. This can be observed 
directly by comparing the visibility curves at each wave¬ 
length (see bottom left panel of Figure 2 ). The shorter 
the wavelength, the faster the visibilities approach 0 as 
the baseline increases, indicating that the emission is get¬ 
ting more and more compact as the wavelength increases. 
Because the central rings are optically thick, this change 
in the derived surface density as a function of wavelength 
could, at least partly, also be due to the difference in opti¬ 
cal depths. However, we explored several models, chang¬ 
ing the dust distribution, but could not find a model that 
reproduces the surface density profile at all wavelengths, 
without invoking a change of the dust distribution as a 
function of the position. 

To build a single model that reproduces the three 
bands simultaneously, we adopt a procedure similar to 
the one used in Pinte et al. (2007) and Gonzalez et al. 
( 2012 ). We use the three previously derived surface den¬ 
sities and assign each of them to a single grain size, which 
emits most at the corresponding wavelength: the sur¬ 
face densities from bands 3, 6 , and 7 were used to de¬ 
fine the distribution of the grains of size 0.4, 0.18, and 
0.11mm respectively. For intermediate sizes, the spatial 
distributions are interpolated in logarithm of the grain 
size. Outside of the range 0.11-0.40 mm, we do not ex¬ 
trapolate the distributions: all the grains smaller than 
0 . 11 mm have the same surface density distribution, in¬ 
dependently of their size. Similarly, all the grains larger 
than 0.40 mm have the same surface density distribu¬ 
tion. Because a given grain size does not emit only in 
one given ALMA band, the final model does not show a 
as good agreement as the individual models only fitting 
one wavelength. However, the emissivity as a function of 
wavelength of a given grain size remains peaked enough 
so that the contamination between bands remains lim¬ 
ited and the general agreement is good at all wavelengths 
(Figure 2 ). The size distribution at any point is no longer 
a power-law. In bottom left panel of Figure 3, we show 
the slope that fits best the local (integrated vertically) 
grain size distribution. The grain distribution has a local 
slope close to -3.5 up to 75 AU, while it approaches -4.5 
between 75 and 120 AU, supporting the idea of depletion 


of large grains. The total dust mass of the best model is 
51O-^M0. 

Table 1 summarizes the properties of the various rings. 
The ring masses range from 7 to 120 M 0 . We use the 
same nomenclature as in ALMA 15 to name the gaps and 
rings. The ring masses are calculated by integrating the 
mass between two successive gaps, and the temperatures 
are the mass averaged dust temperatures. Uncertain¬ 
ties on the gap masses were computed by varying the 
gap depth, i.e. by multiplying the surface density in the 
gap by a constant without varying its shape. This proce¬ 
dure was repeated for various values of the constant until 
the gap was infinitely deep (the surface density virtually 
reaches 0 in the gap) and until it disappeared (the surface 
density in the gap reaches the interpolation between the 
rings). For each of these models, synthetic MCFOST 
images were calculated, visibilities computed using the 
CASA simulator and the total calculated. Due to 
the computing cost of the procedure (dominated by the 
calculations of the visibilities in CASA), a full Bayesian 
analysis taking into account the potential correlations be¬ 
tween the densities in the various gaps and rings is out 
of reach with current computing facilities. Instead, we 
explored the mass of each gap independently. This ap¬ 
proach allowed us to get an estimate of the uncertainties 
on the gaps (and rings) masses. In Table 1, we report 
the interval range where the remains smaller than 
1.5 the minimum y^ value. We find that this threshold 
corresponds to visible differences in the binned visibility 
curves as well as residuals above 3 a in the CLEAN maps. 
Note that the central peak and ring B 1 are optically thick 
even at 3 mm, so the corresponding mass estimates are 
lower limits. 

The dust opacity and spectral index for each ring are 
computed for the integrated dust population in each ring 
(Table 1). The absolute values of the opacity result from 
our choice of dust properties and should be considered 
with an uncertainty of a factor a few. The outer rings dis¬ 
play smaller opacities at 1 mm and spectral index (from 

1 in the central rings to 1.3-1.4 in the two outer rings 
B 6 and B7) due to the steeper grain size distribution 
and relative lack of millimeter-sized dust grains. Note 
that the change in opacity index between the inner and 
outer rings is smaller than the observed change in spec¬ 
tral index (bottom panel of Figure 1), as the low observed 
spectral index in the central regions of the disk results 
from a combination of the lower opacity index and high 
optical depth effects. 

So far in the paper, the iterations on the best model 
have been carried out in the image plane. ALMA pro¬ 
vides images with good fidelity; nevertheless, we also 
check in the uv plane that the final model visibilities 
are consistent with the data. We compare the visibility 
functions of the data sets in all three bands and in the 
direction in u^v space corresponding to the disk major 
and minor axes. To improve SNR, this includes all data 
within ±10° of these angles. The lower panels in Figure 2 
illustrate these functions at each wavelength (left panel) 
and along the minor and major disk axis (right panel), 
showing that there is good agreement between the model 
and observed visibilities. In Figure 2 (top panel), we 
compare the full ALMA image combining bands 6 and 
7 with the model as ’’observed” by the CASA Simula- 
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Figure 3. Structure of the disk model. In all panels, band 3 is shown as red, 6 as green, and 7 as blue lines. Top left: dust surface density 
profiles obtained from inverting the brightness profiles along the disk major axis. The spatial resolution is 10, 4, and 3AU respectively. 
The dashed line represents a tapered-edge model with Rc = 80 AU and 7 = — 0.2 as derived by Kwon et al. (2011, 2015) from CARMA 
1.3 mm data. The gaps and ring identified by ALMA15 are indicated. Top right: vertical optical depth as a function of radius at the three 
wavelengths. Bottom left: exponent of the grain size distribution within a column, using the combined data. The exponent is obtained by 
fitting a power-law to the actual local grain size distribution. Bottom right: map of the gas-to-dust ratio. Only the disk is plotted here, 
not the atmosphere. 


Table 1 

Ring properties. 



Dust 

Mass 

Optically Dust Opaci- 

Dust 

Ring^ 

Mass 

Average 

Tdust 

(K) 

Thick^^ ty*^ (cm^.g“^ 

at of dust) at 

A=3mm A=lmm 

opacity 
Index ^ 

Peak 

>2.3 

630 

X 6.0 

0.9 

B1 

>47 

48 

X 5.7 

1.0 

B2 


36 

marginally 6.6 

1.0 

B3 

soils 

32 

marginally 6.9 

1.1 

B4 

621^0 

25 

6.8 

1.0 

B5 

7+1.7 

-1.5 

22 

7.1 

0.8 

B 6 

lOll,® 

20 

marginally 4.6 

1.4 

B7 

123tlt 

16 

4.9 

1.3 


Notes. 

^ We use the same nomenclature as in ALMA15 to name the rings 
^ Ring masses are calculated by integrating the surface density between 
two successive gaps. Uncertainties correspond to the range of masses 
where ^ 1-5 X miny^ (see the text) 

A ring is considered optically thick if its vertical optical depth is 
larger than 3, and marginally optically thick if it is larger than 1 
^ The dust opacity and opacity index are computed from the integrated 
dust distribution within a ring. 

tor, assuming the uv coverage of the actual observations. 
Much of the faint non-axisymmetric structure seen in 
the image is reproduced by this simulation, indicating 
that this is not inherent to the source. Residuals (top 
right panel) are found to be < 3cr or <1% of the local 
flux/beam in this image. 

4.2. Gap brightness profiles 

The gaps show very similar profiles at bands 6 and 7 
(Figure 3) indicating that they are well resolved at these 


two wavelengths. As a consequence, we used the pro¬ 
file derived from the combined band 6+7 image (with 
the highest signal-to-noise), to measure the properties of 
each gap (Table 2). We define the depth of a gap as the 
ratio of the surface density in the gap and the surface 
density just outside of it (averaged on both sides). The 
gap width is defined as the width at half-depth. These 
values, together with the properties of the bright rings 
(Table 1), can be used to set up hydro dynamical simu¬ 
lations aiming at better understanding the origin of the 
gaps in the disk of HL Tau, although this is outside the 
scope of the present paper. 

The gaps, however, will be ’’filled in” if they are not 
fully spatially resolved, and so this table gives only a 
lower limit to their depth. If we assume that the disk 
surface density was a smooth function before the gaps 
formed, then an estimate of the ’missing mass’ in the 
gaps can be made. The gaps are consistent with a few 
tens of M 0 of solids missing (see Table 2), which is typ¬ 
ical of the core size required for runaway gas accretion, 
leading to gas giants. The uncertainties on the mass of 
the gaps were estimated using a procedure similar to the 
one described in Section 4.1. The density of each rings 
was varied independently until the reached a value 
larger than 1.5 times the minimum y^. 

4.3. Degree of dust settling 

Dust settling in the disk of HL Tau was suggested by 
Brittain et al. (2005), based on the high gas-to-dust ratio 
in the upper layers of the disk, and by Kwon et al. (2011) 
based on 1.3 and 2.6 mm CARMA data and SED, but 
no quantitative estimate was possible due to the limited 
spatial resolution. Because the disk is inclined by ~ 47° 
(with a variation of ±1° between rings, ALMA15), the 
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Table 2 

Gap properties. 


Gap^ 

Position 

(AU) 

Midplane 
Tdust [K] 

Missing 
Dust Mass^ 
(M®) 

Gap 

Depth^ 

Gap 

Width 

(AU) 

DI 

13.2 

45 

to 

•| + 
P ^ 

18 

12 

D2 

32.3 

36 

r\r\ + 6.i 

^^- 4.2 

16 

11 

D3 

42.0 

29 

20ll? 

6.9 

6.6 

D4 

50.0 

24 

7 + 2.4 
' - 5.0 

3.8 

4.5 

D5 

64.2 

23 

q+0.9 

°- 4.7 

8.0 

12 

D6 

73.7 

23 

97+iO 

^‘-17 

12 

8.1 

D7 

91.0 

21 

0(^+ll 

'^'^-20 

11 

9.9 





16 

23 


Notes. 

^ We use the same nomenclature as in ALMA15 to name the gaps 
^ The missing mass is calculated by integrating the difference between 
the surface density in the gap and the interpolated surface density 
between the surrounding two rings. Uncertainties correspond to the 
range of masses where ^ 1-5 X min (see the text) 

^ The gap depth is the ratio of the density in the surrounding rings 
(average of both sides) divided by the density at the minimum of the 
gap. 

° The values of D5+D6 correspond to the assumption that the two 
gaps, D5 and D6, are a large unique gap. 

ALMA images enable us to probe not only the radial 
but also the vertical extension of the various rings. The 
fact that they remain sharp and well-defined at all az¬ 
imuthal angles shows that the dust grains responsible for 
the emission are located in a thin layer in the midplane. 
If the disk were thick, the gap contrast would be reduced 
- particularly along the minor axis of the disk - because 
of projection effects: material in the nearby rings, lo¬ 
cated at significant heights above the midplane, would 
hide the gaps. Figure 4 shows indeed that the appear¬ 
ance of the disk is strongly dependent on the degree of 
dust settling. The variation of the width of the gaps as a 
function of azimuth is directly related to the thickness of 
the emitting region (as well as to the beam shape, which 
is only slightly misaligned with the semi-minor axis). 

A small value of the millimeter dust scale height, himm 
= 0.70 AU, equivalent to cf = 3 x 10“^ produces the 
best agreement with most of the observed features in the 
ALMA image, with rings that are well separated at all 
azimuthal angles. By contrast, in the absence of dust 
settling and for larger values of the a parameter, e.g. 
3 X 10“^, dust of all sizes will be better mixed with the 
gas in a geometrically thicker disk and the millimeter gap 
contrast will be reduced. 

The narrow ring (B5) inside the widest gap (between 
the gaps D5 and D6, ALMA15) shows clear brighten¬ 
ing at the ansae, due to its low optical depth. However, 
the himm = 0.70 AU model is not completely consistent 
with the ALMA observations: this ring seems to disap¬ 
pear along the minor axis (this can also be seen in the 
difference between model and data in the visibility plots 
along the minor axis in Figure 2). This is too large to 
be accounted for by the difference in ring inclination an¬ 
gles. The model with himm = 2.15 AU reproduces this 
behavior but does not provide a strong enough contrast 
in the other rings. This suggests that the degree of ver¬ 
tical stirring of the dust grains may vary as a function of 
the position tof he disk, either due to radial variation of a 
or to local changes in the gas/dust ratio in the midplane 


(which were assumed to be constant with radius here), 
or other stirring mechanism, like a big body in that gap. 

5. THE GAS DISK: AND HGO+. 

HL Tan lies in a relatively dense part of the Taurus 
molecular cloud. Emission from the extended ambient 
cloud contaminates the line profiles near the cloud ve¬ 
locity of 6.8km.s“^ and is resolved out on long baselines. 
However, as noted previously (ALMA15) a spatially com¬ 
pact component is seen in the ALMA images of HCO+ at 
velocities outside that of the narrow ambient emission, 
with a velocity gradient consistent with disk rotation. 
Using the calibrated ALMA Science Verification data 
set covering ^^CO, we have employed CASA CLEAN 
to make a spectral image with a resolution of 0.12 arcsec. 
Like HCO+, this also shows a compact structure in the 
high-velocity ^^CO J=l-0 emission; Eigure 5 compares 
the integrated high-velocity blue- and red-shifted gas in 
the two molecules with the dust emission contours. Both 
lines show relatively bright compact emission coincident 
with the disk, with the highest velocity gradient toward 
the star and oriented along the continuum major axis (or¬ 
thogonal to the outflow direction). The CO image shows 
faint extended red-shifted gas from the outflow to the 
southwest; however, most of the compact high-velocity 
emission appears to be associated with the disk. 

In Eigure 6, we compare position-velocity diagrams 
along the disk major axis in the two molecules, with 
curves representing Keplerian rotation in a geometrically 
thin disk. We use a disk inclination of 47°, and show 
lines for enclosed masses of 1.7 Mq (solid) and 1.0 Mq 
(dashed). Although not formal fits, the high gas veloc¬ 
ities and shape of these PV diagrams suggest that the 
stellar mass is closer to ^1.7M0; this compares with 
the estimates of 1.2 Mq based on HR tracks (Giidel et al. 
2007) and 1.3 Mq based on the approximate extent and 
width of the HCO+ emission line in ALMA15. The gas 
inner radius of 10 AU is an upper limit set by the CO 
surface brightness sensitivity at the highest velocities. 
The outer gas radius is likely a lower limit because of 
the ambient cloud confusion at the low relative veloc¬ 
ities. The location of the second millimeter-grain gap 
(D2) is also indicated in the PV diagrams, and there is 
marginal evidence of an increase in the CO and HCO+ 
surface brightness beyond this radius. Eurther studies 
and full visibility modeling of the gas disk would be pos¬ 
sible with a higher SNR data set at similar resolutions to 
this data, preferably using a molecular transition show¬ 
ing more contrast against the ambient gas. 

6 . DISGUSSION - A GOHERENT MODEL FOR THE HL TAU 

DISK? 

We have produced a radiative transfer model that re¬ 
produces both the multi-wavelength ALMA observations 
and SED of the HL Tan disk (see Appendix A). By 
inverting the brightness maps obtained by ALMA, we 
can reconstruct the millimeter-grain surface density pro¬ 
file of the disk. The central rings (< 30 AU) appear 
marginally optically thick even at 2.9 mm, while the out¬ 
ermost ring becomes optically thin at the longest wave¬ 
length. The comparison of the three ALMA data sets 
between 870 jamm and 3 mm shows that this outer ring 
is also depleted in the larger-sized grains compared to 
the more inner regions, suggesting that grain growth is 
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Figure 4. Effect of dust settling on the appearance of the images at 1 mm. This shows the enhanced contrast as the dust disk gets flatter. 
Upper row: Band 6+7 CLEAN maps, left panel: observations, and three panels on the right: from left to right, model maps without 
settling, with ol = 10“^, and 3 10““^. The corresponding scale height of the grains of 1 mm in size is also indicated on each panel. Middle 

row: cuts of the observed (red) and model (blue) CLEAN maps along the minor axis of the disk (green dashed line in top left panel). 
Bottom row: real part of the Band 7 visibilities (first spectral window) along the disk minor axis, on an expanded scale (red lines are data, 
blue are the models). The visibilities at smaller baselines are not plotted here because they are practically not affected by dust settling. 
The indicated were computed using the full data sets, i.e. all the n, v points in all the spectral windows of the three ALMA bands. 


more efficient in the central regions or that some radial 
migration has already occurred within the disk itself. 

Because our model requires high optical depths at mil¬ 
limeter wavelengths, the calculated midplane tempera¬ 
tures in the gaps (Table 2) are much lower than the tem¬ 
peratures used by Zhang et al. (2015), which assumed 
optically thin emission. As a result, in our model, the 
gaps location does not correspond to the condensation 
fronts of molecules. 

To produce the observed gap/ring contrast, a decrease 
of the density of millimeter-sized grains by a factor of 
>10 is required in the gaps compared to the adjacent 
rings. The corresponding decrement in gas may be much 
more modest and detailed hydro dynamical simulations 
of gas+dust systems are required to determine the mass 
of any planet responsible for such gaps. Numerical pre¬ 


dictions suggest that a planetary core between 30 M 0 
(Paardekooper & Mellema 2004) and 150 M 0 (Fouchet 
et al. 2010) could produce such a decrement in the dust. 
This is similar to the ’missing dust mass’ in these gaps, as 
derived above. The derived masses are also in agreement 
with the maximum planet mass found by Tamayo et al. 
(2015) in their nominal model: 40 M 0 if the planets are 
not in resonance, and 100 M 0 if the outer three planets 
are in a chain of 4:3 resonances. Dedicated HL Tan sim¬ 
ulations by Dipierro et al. (2015) indicate planet masses 
between 60 and 175 M 0 , which is, depending on the ring, 
of the same order or up to a factor ^10 larger than our 
mass estimate. The ’missing dust mass’ may then have 
formed the rocky cores of these planets. 

The well-defined gaps observed by ALMA at all az¬ 
imuthal angles clearly show that the millimeter-emitting 
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Figure 5. Distribution of red and blue-shifted gas traced by (left) J=l-0 and (right) HCO+ J=l-0, compared with the dust 

continuum (shown by contours). The CO image shows integrated emission from -1 —h5km.s“^ and +10 —hl5km.s“^ as cyan and red 
colors respectively. The HCO+ image shows integrated emission from +2 —h5km.s“^ and +9 —hl3km.s“^ as cyan and red colors. The 
contours illustrate the continuum emission in band 3 at levels of 0.08, 0.24, 1.0, and 2.0mjy.beam“^, which correspond approximately 
to the outer edges of the major dust rings B7, B4, Bl, and the bright central region. The HCO+ and continuum data sets are from the 
released ALMA science verification image cubes, and the CO is from the ALMA visibility data, cleaned with a taper of 1600 kA to give a 
resolution of ~0.12arcsec. Axes scale is given in AU and the resolutions of the molecular line images are shown in the lower left corner of 
each panel. 



HLTau HCO+ 



Distance along midplane (AU) 


Figure 6. Position-velocity diagrams of (left) ^^CO J=l-0 and (right) HCO+ J=l-0, taken along the major axis of the disk. The spectra 
± O.lSarcsec on either side of this major axis have been averaged to form these PV diagrams. The velocity extent where ambient cloud 
emission is likely to contaminate the emission is indicated by the dashed white lines (note that the ambient HCO+ line is narrower, so the 
contaminated velocity range is less than for ^^CO). The black lines show the Keplerian velocity for 1.7M© (solid) and 1.0M© (dashed), 
with inner and outer radii of 10-120 AU. On the 1.7 Mq line, the thick section represents the dark gap D2 in the millimeter dust disk (at 
radius of ~35AU). Spectral/spatial resolution is shown upper left. The data have been spectrally smoothed to a resolution of lkm.s“^. 


dust grains are located in a geometrically thin disk, with 
a scale height himm at 100 AU of no more than 2 AU. This 
is at least 5 times smaller than the gas scale height as¬ 
suming hydrostatic equilibrium (see the lower right panel 
in Figure 3). By using a simple settling model, we have 
reproduced the observations with a turbulent viscosity 
coefficient a of 3 x 10“^. Similar values have been de¬ 
rived by Mulders & Dominik (2012) by fitting the median 
SEDs of three samples of protoplanetary disks surround¬ 
ing Herbig stars, T Tauri stars, and brown dwarfs. For 
HL Tau, we obtained this value by direct comparison 
with images. 

The accretion rate of the disk surrounding HL Tau 


was estimated to 8.710“^ MQ.yr”^ from the Bry line lu¬ 
minosity (Beck et al. 2010). Note that the Bry emission 
around HL Tau is extended, so this value should be con¬ 
sidered as an upper limit. The viscosity coefficient can 
be estimated via 

“ 37rc,(r)/i(r)I](r) 

At r = 100 AU, we obtain a ~ 10“^, i.e. at least one or¬ 
der of magnitude larger than the value we estimated from 
the observed physical height of the millimeter grains. 
Both estimates of a depend linearly on the dust-to-gas 
mass ratio, which cannot account for the difference. 
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It is important to note that the two estimates of a do 
not measure the same physical mechanism. The dust disk 
scale height depends on the diffusion coefficient, i.e. on 
the fluctuations of velocities and their correlation time 
(Fromang & Papaloizou 2006), which are then converted 
to an a value assuming a Schmidt number Sc, here set 
to 1.5, while the accretion rate gives an estimate of the 
transport of angular momentum. A larger Sc number, 
i.e. a larger ratio of the momentum and mass diffusivi- 
ties, can explain part of the disagreement. Larger values 
than our assumed Sc = 1.5 have been reported from nu¬ 
merical simulations e.g. 2.8 in Fromang & Papaloizou 
(2006) and 11 in Carballido et al. (2005). A simple a- 
prescription for the viscosity might also not be accurate 
enough to describe both mechanisms simultaneously. In 
particular, the a coefficient is very likely not constant 
throughout the disk and the estimation of M might be 
linked to different parts of the disk. The dust diffusion 
process was also found to be highly anisotropic, with 
larger values of Sc in the radial direction than in the 
vertical direction (e.g. Johansen et al. 2006b). 

The SED of HL Tau - particularly the far-infrared 
fluxes and the observed extinction to the star - cannot be 
reproduced with a hydrostatically supported disk. Pos¬ 
sible explanations include a large envelope, or alterna¬ 
tively, a puffed-up disk atmosphere. The steep redden¬ 
ing gradient within 1 arcsec of the disk noted by Close 
et al. (1997) suggests that most of the extinction is in 
the disk rather than an extended envelope. As the gas 
and micron-sized grains are closely coupled, if the gas 
scale height at 100 AU was 30 AU, it would be possible 
to fit the SED and get such a large extinction gradient. 
This may be infalling gas, or part of a photoevaporative 
or MHD wind (e.g. Gorti et al. 2015). 

With the assumptions of our model, the gas-to-dust 
mass ratio in the disk midplane is « 6 at a radius of 
100 AU, if we integrate the dust size distribution up to 
3 mm. If as expected, grains have grown further in the 
midplane, then the actual gas-to-dust ratio may be even 
lower. As an indication, and assuming our model can 
be extrapolated to larger sizes, gas-to-dust ratio of 1 at 
100 AU is reached for Umax = 10 cm. 

Because the millimeter dust disk is also very geomet¬ 
rically thin, these results provide direct constraints on 
the models of solid growth in the disk midplane (Youdin 
& Shu 2002; Johansen et al. 2006a). The increased dust 
density leads to a vertical rotational velocity gradient 
near the midplane because of the feedback (collision) of 
dust onto the gas. The shear on the gas rotational ve¬ 
locity increases as the dust settles and the gas/dust ratio 
decreases, leading eventually to Kelvin-Helmholtz insta¬ 
bilities, which in turn will prevent the dust grains from 
settling further. It has been suggested that, in this con¬ 
figuration, streaming instabilities will lead to very rapid 
clumping of dust particles (Youdin & Goodman 2005; 
Johansen & Youdin 2007), possibly providing a mecha¬ 
nism to grow boulders rapidly. The increased gas/dust 
ratio and flat dust sub-disk can also lead to ’’pebble ac¬ 
cretion” (Lambrechts & Johansen 2012), significantly re¬ 
ducing the growth timescale of a planetary core massive 
enough to trigger runaway gas accretion. In the calcula¬ 
tions, turbulent diffusion and shear instabilities limit the 
amount of dust accumulation in the midplane, therefore 


limiting the growth of the big bodies. Our results, al¬ 
though limited to particle sizes of a few millimeters, give 
specific indications regarding the size of the dust sub-disk 
and the decrease in the gas/dust ratio. 

7. CONCLUSIONS 

We have produced a model that reproduces the ALMA 
high-resolution observations of HL Tau at the three wave¬ 
lengths. Detailed features including the rings, the gaps, 
their contrast and emissivity index can be modeled. The 
results indicate that the gaps are devoid of dust by at 
least a factor of 10 compared to the adjacent rings. The 
dust masses in the rings range from a few to ^ 100 M 0 , 
and we estimate that if the disk originally had a smooth 
radial distribution, up to 40 of dust has been re¬ 
moved to create each of the two largest gaps. 

Our modeling shows that the increased optical depth in 
the central rings cannot explain alone the change in ob¬ 
served spectral index with radius. A change in the dust 
properties is also required, which indicates that larger 
grains are present closer to the star. This suggests that 
grain growth is faster in the central regions and/or that 
radial migration is occurring in the disk. The high opti¬ 
cal depth of the disk in our model also means that the 
temperature in the gaps is significantly lower than what 
was assumed by Zhang et al. (2015) and that the gaps do 
not seem to be coincident with molecular condensation 
fronts. 

Interestingly, the sharp rings observed at all azimuthal 
angles clearly show that the disk emitting at millime¬ 
ter wavelengths is geometrically thin: of order 1 AU at 
100 AU. We interpret this as a clear sign of dust verti¬ 
cal settling. Assuming a standard dust settling model, 
this implies that the coefficient of turbulent viscosity is 
of the order a few 10“^, between one and two orders of 
magnitude smaller than what is required to account for 
the estimated accretion on the protostar. Because most 
dust mass is in the larger grains, this would indicate that 
the gas/dust ratio is about 5 in these rings. Eor the first 
time (to our knowledge), we provide numbers obtained 
directly from observations for the thickness and the de¬ 
crease in the gas/dust ratio in the disk mid-plane. By 
contrast, in order to fit the SED, the disk thickness in 
small grains and gas may be significantly larger than hy¬ 
drostatic equilibrium. 

The parameters of the millimeter-grain disk are inde¬ 
pendent of the origins of the gaps and rings, but they 
will no doubt provide significant constraints for planet 
formation models and/or HD or MHD models that will 
aim to reproduce them. 

This paper makes use of the following ALMA data: 
ADS/JAO.ALMA#2011.0.00015.SV. ALMA is a part¬ 
nership of ESO (representing its member states), NSE 
(USA) and NINS (Japan), together with NRG (Ganada) 
and NSG and ASIAA (Taiwan), in cooperation with the 
Republic of Ghile. The Joint ALMA Observatory is op¬ 
erated by ESO, AUI/NRAO, and NAOJ. The National 
Radio Astronomy Observatory is a facility of the Na¬ 
tional Science Eoundation operated under cooperative 
agreement by Associated Universities, Inc. I.d.G. ac¬ 
knowledges support from MIGINN (Spain) AYA2011- 
30228-G03 grant (including EEDER funds). We thank 
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C. Dougados, G. Lesur, and H. Klahr for useful discus¬ 
sions. The research leading to these results has received 
funding from the European Union Seventh Framework 
Programme FP7-2011 under grant agreement no 284405. 


APPENDIX 

A. FITTING SIMULTANEOUSLY THE SED AND 
MILLIMETER MAPS 

Small (sub-micron) sized dust grains will be well-mixed 
with the gas. However, the simple model of a thin hydro¬ 
static disk with a gas (and sub-micron dust) scale height 
ho ~10AU cannot reproduce the SED, particularly in 
the far-infrared. Several authors have noted that an ad¬ 
ditional component is needed to fit HL Tan, and generally 
this has been assumed to be a spherical infalling envelope 
(Men’shchikov et al. 1999; Robitaille et al. 2007; Kwon 
et al. 2011). However, the ALMA data (ALMA15) also 
show that this system has a broad CO bipolar outflow 
with an opening angle of ^ 80° - indicating that such an 
envelope would not be a spherical structure. The incli¬ 
nation of the system suggests that we are looking along 
a line of sight somewhere between this outflow cone and 
the ho = 10 AU gas scale height proposed above. Several 
high-density molecular tracers (CN, HCN, and HCO+) 
show narrow line-of-sight absorption against the bright 
disk continuum, and the line profiles suggest the gas may 
be blue-shifted by ^ 0.5km.s“^ (ALMA15). This would 
imply that the gas high up in the disk atmosphere is 
actually outflowing rather than infalling. The line-of- 
sight gas column to the star has also been probed by 
infrared CO lines (Brittain et al. 2005), and a compari¬ 
son of the derived gas and dust column density was used 
to show that this material is dust-depleted, i.e. has a 
gas/dust ratio of >100. Because the photometry from 
the optical to mid-infrared regime suffers from high opti¬ 
cal depth extinction, it is not possible to assess the exact 
nature of the dust structure responsible for this extinc¬ 
tion. Several dust spatial distributions can account for 
the optical and near-infrared extinction and mid-infrared 
excess above the disk emission. We chose to add a thick 
disk atmosphere composed of interstellar-like dust grains 
(same composition as in the disk but with Umax = 1 /im). 
We find a good agreement with the observed SED for 
ho = 30 AU, a surface density slope of -0.75 and dust 
mass of 2 . 81 O“^M 0 . The resulting model SED is pre¬ 
sented in Figure 7. It does not require the addition of 
a spherical envelope, but does require a disk in the sub¬ 
micron sized grains that is considerably thicker than hy¬ 
drostatic equilibrium. Because of its low mass, compared 
to the disk, and small opacity at millimeter wavelength, 
this extra component is undetectable in the synthetic 
ALMA maps and does not affect the results presented in 
this paper. 

B. HYDROSTATIC SCALE HEIGHT 

In the modeling presented in this paper, we have as¬ 
sumed a fixed gas vertical structure. Figure 8 presents 
the calculated midplane temperatures compared to the 
temperature corresponding to the gas scale height of the 
best model. The agreement is very good almost every¬ 
where in the disk, except at the very inner edge where 



Figure 7. Spectral energy distribution of HL Tan. Photometry 
data points are from Robitaille et al. (2007), Men’shchikov et al. 
(1999), and ALMA15. The triangles show the fluxes obtained by 
ALMA. Red solid line: disk + atmosphere, dashed blue line: disk 
only. 



Figure 8. Gas vertical scale height. The full red line represents 
the hydrostatic scale height computed from the midplane temper¬ 
ature of the best model. The dashed black line represents the gas 
scale height used as an input for the model, i.e. 10 AU at 100 AU. 
The shaded area represents variations of the scale height of 10 %, 
i.e. between 9 and 11 AU at 100 AU. 

stellar radiation heats the midplane directly. In the outer 
disk, the actual disk temperature oscillates due to the al¬ 
ternating rings and gaps, but the variation remains lim¬ 
ited and corresponds to a change in the scale height of 
about 5%. 
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